Here we fit a survival model to lifetimes of patients with Acute Myeloid Leukemia.
data(Leuk)
For survival models we need to scale the time for numerical stability.
## Set survival time as year
Leuk$time <- Leuk$time / 365.25
Let’s look at some descriptive statistics.
## Data summary
round(sapply(Leuk[, c(1, 2, 5:8)], summary), 2)
## time cens age sex wbc tpi
## Min. 0.00 0.00 14.00 0.00 0.00 -6.09
## 1st Qu. 0.11 1.00 49.00 0.00 1.80 -2.70
## Median 0.51 1.00 65.00 1.00 7.90 -0.37
## Mean 1.46 0.84 60.73 0.52 38.59 0.34
## 3rd Qu. 1.47 1.00 74.00 1.00 38.65 2.93
## Max. 13.63 1.00 92.00 1.00 500.00 9.55
## a visualization of the data locations and boundary
plot(nwEngland)
with(Leuk,
points(xcoord, ycoord,
pch = 16,
col = cens + 1,
cex = 0.1 + log(time/4)))
## plot the Kaplan-Mayer estimator
km <- survfit(Surv(time, cens) ~ sex, Leuk)
par(mar = c(2.5, 2.5, 0.5, 0.5),
mgp = c(1.5, 0.5, 0), las = 1)
plot(km, conf.int = TRUE,
col = 2:1, bty = "n")
legend('topright', c('female', 'male'),
lty = 1, col = 2:1,
bty = "n")
Now we build the models. We will consider a non-spatial model and a
Matern field model.
## model formulae without u()
form0 <- inla.surv(time, cens) ~
sex + age + wbc + tpi
## fitting the model without u()
fit0 <- inla(
formula = form0,
family = "weibullsurv",
data = Leuk)
fit0$summary.fixed
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -2.038354804 0.1422330553 -2.317285421 -2.038344445 -1.75948333
## sex 0.067865055 0.0676957061 -0.064879285 0.067864975 0.20060985
## age 0.030148635 0.0020718871 0.026086128 0.030148556 0.03421160
## wbc 0.002913092 0.0004532199 0.002024694 0.002912978 0.00380214
## tpi 0.025121025 0.0089954746 0.007481389 0.025121164 0.04275987
## mode kld
## (Intercept) -2.038344434 6.806084e-11
## sex 0.067864975 5.526746e-11
## age 0.030148556 6.167659e-11
## wbc 0.002912977 6.865013e-11
## tpi 0.025121165 5.512512e-11
fit0$summary.hyperpar
## mean sd 0.025quant 0.5quant
## alpha parameter for weibullsurv 0.5791636 0.01488211 0.550262 0.5790751
## 0.975quant mode
## alpha parameter for weibullsurv 0.6085885 0.5788969
Now we set up the mesh for the Matern field.
## 1: build a mesh: for u(s_0), s_0 are mesh nodes
## not a very good one for this problem:
mesh0 <- fm_mesh_2d(
nwEngland, ## where
max.edge = c(0.05, 0.2) ## (max) triangle sizes
)
## show it
plot(mesh0, asp = 1)
plot(nwEngland, add = TRUE, lwd = 2)
## better: non-convex boundaries around
bnd1 <- fm_nonconvex_hull(
x = nwEngland,
convex = 0.03,
concave = 0.1,
resol = 50)
bnd2 <- fm_nonconvex_hull(
x = nwEngland,
convex = 0.25)
## build a mesh accounting for boundaries
mesh <- fm_mesh_2d(
boundary = list(bnd1, bnd2),
max.edge = c(0.05, 0.2),
cutoff = 0.02)
## visualize it
plot(mesh)
plot(nwEngland, add = TRUE, lwd = 2)
## 2. build a projector matrix
data.proj <- fm_evaluator(
mesh = mesh,
loc = cbind(Leuk$xcoord,
Leuk$ycoord)
)
str(data.proj)
## List of 6
## $ x : NULL
## $ y : NULL
## $ lattice: NULL
## $ loc : num [1:1043, 1:2] 0.205 0.286 0.176 0.245 0.327 ...
## $ proj :List of 4
## ..$ t : int [1:1043, 1] 1455 615 303 1394 1138 1064 669 520 519 814 ...
## ..$ bary: num [1:1043, 1:3] 0.2726 0.4114 0.0078 0.0296 0.3979 ...
## ..$ A :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:3129] 194 60 271 442 663 682 1020 346 708 195 ...
## .. .. ..@ p : int [1:893] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. ..@ Dim : int [1:2] 1043 892
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:3129] 0.03384 0.08384 0.00929 0.09469 0.09432 ...
## .. .. ..@ factors : list()
## ..$ ok : logi [1:1043] TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ crs : NULL
## - attr(*, "class")= chr "fm_evaluator"
## 3. build the SPDE model object
spde <- inla.spde2.pcmatern(
mesh = mesh, ## where to build on
prior.range = c(0.05, 0.01), ## P(range < 0.05) = 0.01
prior.sigma = c(1, 0.01)) ## P(sigma > 1) = 0.01
## model formulae adding u(): add f() model term
## index needed for f(): NA's as we will use A.local
Leuk$spatial <- rep(NA, nrow(Leuk))
## projector matrix (from mesh nodes to obs. locations)
form1 <- update(
form0,
. ~ . + f(spatial, model = spde,
A.local = data.proj$proj$A))
## mode fit with u()
fit1 <- inla(
formula = form1,
family = "weibullsurv",
data = Leuk
)
fit1$summary.fix
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -2.170303113 0.2056588067 -2.569461854 -2.173013319 -1.753059507
## sex 0.072267360 0.0691664958 -0.063405229 0.072282232 0.207855211
## age 0.032965443 0.0022352896 0.028594357 0.032961435 0.037359091
## wbc 0.003074293 0.0004579608 0.002176741 0.003074126 0.003972794
## tpi 0.024619020 0.0098421183 0.005347963 0.024607647 0.043954980
## mode kld
## (Intercept) -2.172411962 3.416398e-08
## sex 0.072282287 5.873656e-11
## age 0.032961341 1.059823e-09
## wbc 0.003074125 8.076372e-11
## tpi 0.024607640 3.542345e-10
fit1$summary.hyperpar
## mean sd 0.025quant 0.5quant
## alpha parameter for weibullsurv 0.5993735 0.01600970 0.5684616 0.5991635
## Range for spatial 0.3103588 0.15551112 0.1139983 0.2762600
## Stdev for spatial 0.2925458 0.07319816 0.1741384 0.2839863
## 0.975quant mode
## alpha parameter for weibullsurv 0.6314929 0.5987528
## Range for spatial 0.7087770 0.2203445
## Stdev for spatial 0.4604135 0.2678006
## build a projector for a grid
(bbnw <- bbox(nwEngland))
## min max
## x -0.004 0.832
## y -0.047 1.077
r0 <- diff(range(bbnw[1, ])) / diff(range(bbnw[2, ]))
grid.proj <- fm_evaluator(
mesh = mesh,
xlmim = bbnw[1, ],
ylim = bbnw[2, ],
dims = c(200 * r0, 200)
)
str(grid.proj)
## List of 6
## $ x : num [1:149] -0.253 -0.244 -0.235 -0.226 -0.217 ...
## $ y : num [1:200] -0.047 -0.0414 -0.0357 -0.0301 -0.0244 ...
## $ lattice:List of 6
## ..$ dims: int [1:2] 149 200
## ..$ x : num [1:149] -0.253 -0.244 -0.235 -0.226 -0.217 ...
## ..$ y : num [1:200] -0.047 -0.0414 -0.0357 -0.0301 -0.0244 ...
## ..$ loc : num [1:29800, 1:2] -0.253 -0.244 -0.235 -0.226 -0.217 ...
## ..$ segm:List of 5
## .. ..$ loc : num [1:694, 1:3] -0.253 -0.244 -0.235 -0.226 -0.217 ...
## .. ..$ idx : int [1:694, 1:2] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ grp : int [1:694] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ is.bnd: logi [1:694] TRUE TRUE TRUE TRUE TRUE TRUE ...
## .. ..$ crs : NULL
## .. ..- attr(*, "class")= chr [1:2] "fm_segm" "inla.mesh.segment"
## ..$ crs : NULL
## ..- attr(*, "class")= chr [1:2] "fm_lattice_2d" "inla.mesh.lattice"
## $ loc : NULL
## $ proj :List of 4
## ..$ t : int [1:29800, 1] NA NA NA NA NA NA NA NA NA NA ...
## ..$ bary: num [1:29800, 1:3] NA NA NA NA NA NA NA NA NA NA ...
## ..$ A :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:77295] 11178 11326 11327 11328 11475 11476 11477 11478 11624 11625 ...
## .. .. ..@ p : int [1:893] 0 467 1109 1269 1396 1709 2255 2465 2837 2929 ...
## .. .. ..@ Dim : int [1:2] 29800 892
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:77295] 0.0227 0.1002 0.0595 0.0188 0.137 ...
## .. .. ..@ factors : list()
## ..$ ok : logi [1:29800] FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ crs : NULL
## - attr(*, "class")= chr "fm_evaluator"
## evaluate the field
spat.m <- fm_evaluate(
grid.proj,
field = fit1$summary.random$spatial$mean)
str(spat.m)
## num [1:149, 1:200] NA NA NA NA NA NA NA NA NA NA ...
image(grid.proj$x,
grid.proj$y,
spat.m,
asp = 1)
plot(nwEngland, add = TRUE, border = "lightgrey")
## improve the plot: set as NA outside the map
ov <- over(SpatialPoints(grid.proj$lattice$loc),
nwEngland)
spat.m[is.na(ov)] <- NA
image(grid.proj$x,
grid.proj$y,
spat.m,
asp = 1)
plot(nwEngland, add = TRUE)
## improve the plot
i.ok <- which(!is.na(spat.m))
xyz <- data.frame(
x = grid.proj$lattice$loc[i.ok, 1],
y = grid.proj$lattice$loc[i.ok, 2],
z = as.vector(spat.m)[i.ok]
)
gg0 <- ggplot() + theme_minimal()
nwEngland_geom <- list()
for (i in 1:length(nwEngland@polygons[[1]]@Polygons)){
nwEngland_geom <- c(nwEngland_geom, list(nwEngland@polygons[[1]]@Polygons[[i]]@coords))
}
gg.m <- gg0 +
geom_contour_filled(
data = xyz,
aes(x = x, y = y, z = z)) +
scale_fill_distiller(super = metR::ScaleDiscretised, palette = "Spectral") +
geom_sf(data = st_as_sf(st_sfc(st_polygon(nwEngland_geom))),
fill = 'transparent', lwd = 0.3, col = "black") +
theme(legend.title = element_blank(),
legend.key.height = unit(1, "null")) +
guides(colour = guide_colourbar(position = "bottom"))
gg.m
xyz$sd <- as.vector(
fm_evaluate(
grid.proj,
field = fit1$summary.random$spatial$sd
))[i.ok]
summary(xyz)
## x y z sd
## Min. :-0.0005665 Min. :-0.04135 Min. :-0.543314 Min. :0.1999
## 1st Qu.: 0.2338393 1st Qu.: 0.26365 1st Qu.:-0.124188 1st Qu.:0.2306
## Median : 0.3871046 Median : 0.46699 Median :-0.006666 Median :0.2505
## Mean : 0.3919586 Mean : 0.46881 Mean :-0.030992 Mean :0.2533
## 3rd Qu.: 0.5493855 3rd Qu.: 0.65903 3rd Qu.: 0.113343 3rd Qu.:0.2766
## Max. : 0.8288693 Max. : 1.07135 Max. : 0.352537 Max. :0.3326
gg.sd <- gg0 +
geom_contour_filled(
data = xyz,
aes(x = x, y = y, z = sd)) +
scale_fill_distiller(super = metR::ScaleDiscretised, palette = "Spectral") +
geom_sf(data = st_as_sf(st_sfc(st_polygon(nwEngland_geom))),
fill = 'transparent', lwd = 0.3, col = "black") +
theme(legend.title = element_blank(),
legend.key.height = unit(1, "null")) +
guides(colour = guide_colourbar(position = "bottom"))
gg.sd
## visualize the spatial effect and uncertainty
grid.arrange(
gg.m,
gg.sd,
nrow = 1
)
## Example 2 - spatio-temporal model for precipitation We analyse daily
rainfall in Parana State. It does not rain everyday so we usually need a
hurdle model.
data(PRprec)
head(PRborder,2)
## Longitude Latitude
## [1,] -54.60785 -25.44626
## [2,] -54.59983 -25.43444
border.ll <- SpatialPolygons(list(Polygons(list(
Polygon(PRborder)), '0')),
proj4string = CRS("+proj=longlat +datum=WGS84"))
border <- spTransform(border.ll,
CRS("+proj=utm +units=km +zone=22 +south"))
bbox(border)
## min max
## x 136.0091 799.8802
## y 7044.8116 7509.5597
apply(bbox(border), 1, diff)
## x y
## 663.8711 464.7481
PRprec[1:3, 1:8]
## Longitude Latitude Altitude d0101 d0102 d0103 d0104 d0105
## 1 -50.8744 -22.8511 365 0 0 0 0.0 0
## 3 -50.7711 -22.9597 344 0 1 0 0.0 0
## 4 -50.6497 -22.9500 904 0 0 0 3.3 0
loc.ll <- SpatialPoints(PRprec[,1:2], border.ll@proj4string)
loc <- spTransform(loc.ll, border@proj4string)
#For illustration we use the first 8 days
m <- 8
days.idx <- 3 + 1:m
z <- as.numeric(PRprec[, days.idx] > 0)
table(z)
## z
## 0 1
## 3153 1719
y <- ifelse(z == 1, unlist(PRprec[, days.idx]), NA)
table(is.na(y))
##
## FALSE TRUE
## 1719 3209
mesh <- inla.mesh.2d(loc, max.edge = 200, cutoff = 35,
offset = 150)
par(mar = c(0,0,0,0))
plot(mesh)
plot(loc, add = TRUE, col = "blue", pch = 20, cex = 0.5)
#plot(loc, add = TRUE, col = "blue", pch = 20, cex = 0.5)
plot(border, col = NA, border = "red", add = TRUE)
prange <- c(100, 0.5)
psigma <- c(0.5, 0.5)
spde <- inla.spde2.pcmatern(
mesh, prior.range = prange, prior.sigma = psigma)
rhoprior <- list(rho = list(prior = 'pc.cor1',
param = c(0.5, 0.7)))
bprior <- list(prior = 'gaussian', param = c(0,1))
pcgprior <- list(prior = 'pc.gamma', param = 1)
n <- nrow(PRprec)
#spatio-temporal coords
stcoords <- kronecker(matrix(1, m, 1), coordinates(loc))
A <- inla.spde.make.A(mesh = mesh, loc = stcoords,
group = rep(1:m, each = n))
dim(A) == (m * c(n, spde$n.spde))
## [1] TRUE TRUE
field.z.idx <- inla.spde.make.index(name = 'x',
n.spde = spde$n.spde, n.group = m)
field.zc.idx <- inla.spde.make.index(name = 'xc',
n.spde = spde$n.spde, n.group = m)
field.y.idx <- inla.spde.make.index(name = 'u',
n.spde = spde$n.spde, n.group = m)
stk.z <- inla.stack(
data = list(Y = cbind(as.vector(z), NA), link = 1),
A = list(A, 1),
effects = list(field.z.idx, z.intercept = rep(1, n * m)),
tag = 'zobs')
stk.y <- inla.stack(
data = list(Y = cbind(NA, as.vector(y)), link = 2),
A = list(A, 1),
effects = list(c(field.zc.idx, field.y.idx),
y.intercept = rep(1, n * m)),
tag = 'yobs')
stk.zp <- inla.stack(
data = list(Y = matrix(NA, ncol(A), 2), link = 1),
effects = list(field.z.idx, z.intercept = rep(1, ncol(A))),
A = list(1, 1),
tag = 'zpred')
stk.yp <- inla.stack(
data = list(Y = matrix(NA, ncol(A), 2), link = 2),
A = list(1, 1),
effects = list(c(field.zc.idx, field.y.idx),
y.intercept = rep(1, ncol(A))),
tag = 'ypred')
stk.all <- inla.stack(stk.z, stk.y, stk.zp, stk.yp)
cff <- list(list(), list(hyper = list(prec = pcgprior)))
cg <- list(model = 'ar1', hyper = rhoprior)
formula.joint <- Y ~ -1 + z.intercept + y.intercept +
f(x, model = spde, group = x.group, control.group = cg) +
f(xc, copy = "x", fixed = FALSE, group = xc.group,
hyper = list(beta = bprior)) +
f(u, model = spde, group = u.group, control.group = cg)
# Initial values of parameters - for illustration
ini.jo <- c(-0.047, 5.34, 0.492, 1.607, 4.6, -0.534, 1.6, 0.198)
res.jo <- inla(formula.joint, family = c("binomial", "gamma"),
data = inla.stack.data(stk.all), control.family = cff,
control.predictor = list(A = inla.stack.A(stk.all),
link = link),
control.compute = list(dic = TRUE, waic = TRUE,
config = TRUE), control.inla = list(int.strategy = "eb"),
control.mode = list(theta = ini.jo, restart = TRUE))
formula.zy <- Y ~ -1 + z.intercept + y.intercept +
f(x, model = spde, group = x.group, control.group = cg) +
f(u, model = spde, group = u.group, control.group = cg)
# Initial values of parameters - for illustration
ini.zy <- c(-0.05, 5.3, 0.5, 1.62, 4.65, -0.51, 1.3)
res.zy <- inla(formula.zy, family = c("binomial", "gamma"),
data = inla.stack.data(stk.all), control.family = cff,
control.predictor = list(A =inla.stack.A(stk.all),
link = link),
control.compute=list(dic = TRUE, waic = TRUE,
config = TRUE), control.inla = list(int.startegy = "eb"),
control.mode = list(theta = ini.zy, restart = TRUE))
## Warning in ctrl_check(x = x, the_type = the_type): Control name(s) 'int.startegy' appears to be invalid for `control.inla` and will be removed.
## Valid names are:
## adapt.hessian.max.trials, adapt.hessian.mode, adapt.hessian.scale,
## adaptive.max, adjust.weights, b.strategy, cmin, compute.initial.values,
## constr.marginal.diagonal, control.vb, cpo.diff, cutoff, diagonal,
## diff.logdens, dz, fast, force.diagonal, global.node.degree,
## global.node.factor, h, hessian.correct.skewness.only, huge,
## improved.simplified.laplace, int.design, int.strategy, interpolator,
## lincomb.derived.correlation.matrix, linear.correction, mode.known,
## npoints, num.gradient, num.hessian, numint.abserr, numint.maxfeval,
## numint.relerr, optimise.strategy, optimiser, parallel.linesearch,
## print.joint.hyper, reordering, restart, skip.configurations, stencil,
## step.factor, step.len, strategy, stupid.search, stupid.search.factor,
## stupid.search.max.iter, tolerance, tolerance.f, tolerance.g,
## tolerance.step, tolerance.x, use.directions, verbose
formula.sh <- Y ~ -1 + z.intercept + y.intercept +
f(x, model = spde, group = x.group, control.group = cg) +
f(xc, copy = "x", fixed = FALSE, group = xc.group)
# Initial values of parameters - for illustration
ini.sh <- c(-0.187, 5.27, 0.47, 1.47, 0.17)
res.sh <- inla(formula.sh, family = c("binomial", "gamma"),
data = inla.stack.data(stk.all), control.family = cff,
control.predictor = list(
A = inla.stack.A(stk.all), link = link),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE,
config = TRUE), control.inla = list(int.strategy = "eb"),
control.mode = list(theta = ini.sh, restart = TRUE))
getfit <- function(r) {
fam <- r$dic$family
data.frame(dic = tapply(r$dic$local.dic, fam, sum),
waic = tapply(r$waic$local.waic, fam, sum))
}
rbind(separate = getfit(res.jo),
joint = getfit(res.zy),
oshare = getfit(res.sh))[c(1, 3, 5, 2, 4, 6),]
## dic waic
## separate.1 5081.742 5073.134
## joint.1 5087.927 5078.164
## oshare.1 5087.088 5078.020
## separate.2 11278.805 11291.715
## joint.2 11300.193 11314.125
## oshare.2 11456.897 11457.523
idx.z <- inla.stack.index(stk.all, 'zobs')$data
idx.y <- inla.stack.index(stk.all, 'yobs')$data
idx.zp <- inla.stack.index(stk.all, 'zpred')$data
idx.yp <- inla.stack.index(stk.all, 'ypred')$data
wh <- apply(bbox(border), 1, diff)
nxy <- round(300 * wh / wh[1])
pgrid <- inla.mesh.projector(mesh, xlim = bbox(border)[1, ],
ylim = bbox(border)[2, ], dims = nxy)
ov <- sp::over(SpatialPoints(pgrid$lattice$loc,
border@proj4string), border)
df = st_as_sf(data.frame(pgrid$lattice$loc), coords=c(1,2), crs=border@proj4string)
ov <- lengths(st_intersects(df, st_as_sf(border))) > 0
id.out <- which(ov == FALSE)
stpred <- matrix(res.jo$summary.fitted.values$mean[idx.zp],
spde$n.spde)
stpred_y <- matrix(res.jo$summary.fitted.values$mean[idx.yp],
spde$n.spde)
par(mfrow = c(1, 1), mar =c(0, 0, 0, 0))
for (j in 1:m) {
pj <- inla.mesh.project(pgrid, field = stpred[, j])
pj[id.out] <- NA
image.plot(x = pgrid$x, y = pgrid$y, z = pj,
zlim = c(0, 1), col = viridis(50), breaks = seq(0,1,length = 51),
main = paste0("Day ", j, " - P(rain)"), ylab = "", xlab = "")
}
par(mfrow = c(1, 1), mar =c(0, 0, 0, 0))
for (j in 1:m) {
pj <- inla.mesh.project(pgrid, field = stpred_y[, j])
pj[id.out] <- NA
image.plot(x = pgrid$x, y = pgrid$y, z = pj,
zlim = c(0, 1), col = viridis(50), breaks = seq(0,30,length = 51),
main = paste0("Day ", j, " - Predicted rainfall"), ylab = "", xlab = "")
}
Taken from https://github.com/inlabru-org/inlabru/blob/devel/vignettes/articles/svc.Rmd.
The model used to analyze these data assumes that counts come from a negative binomial distribution with an expected count and dispersion parameter. The expected count has a log-linear predictor:
\[\log(\lambda_{st}) = \kappa_{s} + \alpha_s + \epsilon_{s} \log[\text{Effort}_{st}] + \tau_{s} \text{Year}_{st}\]
where the natural log of expected count, \(\log(\lambda_{st})\), at site \(s\) during year \(t\), is modeled with a zero-centered, normally distributed intercept per site, \(\kappa_s\), a spatially varying intercept, \(\alpha_s\), a spatially varying effect of the log of count effort in hours, \(\epsilon_s\), and a spatially varying linear effect of year, \(\tau_s\). The spatially structured effects are modeled as Gaussian random fields with Matérn covariance functions with range and variance parameters.
Model parameters \(\kappa_s\), \(\alpha_s\), \(\epsilon_s\) and \(\tau_s\) are analogous to those in Meehan et al. (2019). For example, \(\kappa_s\) is included to account for site-level differences in counts, possibly due to habitat availability or observer experience. \(\alpha_s\) can be interpreted as an effort-corrected abundance index at year zero. \(\epsilon_s\) is the exponent for a power-law effort-correction function. And \(\tau_s\) is the long-term temporal trend at a given site.
We define a coordinate reference system (CRS) for spatial analysis and create a base map for later use. The CRS uses the USA Contiguous Albers Equal-Area Conic projection, and is identified by the EPSG code 6703. We modify the CRS slightly to have units of kilometers, so that distances between widespread count sites are not especially large numbers (Krainski et al. 2018).
# define a crs
epsg6703km <- paste(
"+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5",
"+lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83",
"+units=km +no_defs"
)
# make a base map
states <- maps::map("state", plot = FALSE, fill = TRUE) %>%
sf::st_as_sf() %>%
filter(ID %in% c(
"texas", "oklahoma", "kansas", "missouri",
"arkansas", "louisiana"
)) %>%
sf::st_make_valid() %>%
sf::st_transform(epsg6703km) %>%
sf::st_make_valid()
Next we import some bird count data from the GitHub repository associated with Meehan et al. (2019), and turn the data set into spatially referenced points. We use a subset of the data (30 years, 6 US states) for this analysis to reduce computing time (~ 1 min). Note that site selection and zero filling, important components of trend analyses, have already been conducted and this is the resulting data set.
The inlabru package contains a pregenerated version of
the data subset, called robins_subset, that can be accessed
with data(robins_subset), to avoid accessing the full data
online. The following code was used to generate the subset:
robins_subset <- read.csv(paste0(
"https://raw.github.com/tmeeha/inlaSVCBC",
"/master/code/modeling_data.csv"
)) %>%
select(
circle, bcr, state, year, std_yr, count, log_hrs,
lon, lat, obs
) %>%
mutate(year = year + 1899) %>%
filter(
state %in% c(
"TEXAS", "OKLAHOMA", "KANSAS", "MISSOURI",
"ARKANSAS", "LOUISIANA"
),
year >= 1987
)
We load the data, filter out observation sites with less than 20
years of data, add index variables to uniquely index site and year
information, and transform the coordinates to the
epsg6703km CRS:
data(robins_subset)
count_dat <- robins_subset %>%
mutate(site_idx = as.numeric(factor(paste(circle, lon, lat)))) %>%
group_by(site_idx) %>%
mutate(n_years = n()) %>%
filter(n_years >= 20) %>%
ungroup() %>%
mutate(
std_yr = year - max(year),
obs = seq_len(nrow(.)),
site_idx = as.numeric(factor(paste(circle, lon, lat))),
year_idx = as.numeric(factor(year)),
site_year_idx = as.numeric(factor(paste(circle, lon, lat, year)))
) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>%
st_transform(epsg6703km) %>%
mutate(
easting = st_coordinates(.)[, 1],
northing = st_coordinates(.)[, 2]
) %>%
arrange(circle, year)
A rough view of changes in robin relative abundance, which does not account for variation in count effort, can be seen by plotting raw counts per site and year.
# map it
ggplot() +
geom_sf(
data = count_dat %>% filter(year_idx %in% seq(1, 30, 3)),
aes(col = log(count + 1))
) +
geom_sf(data = states, fill = NA) +
coord_sf(datum = NA) +
facet_wrap(~year) +
scale_color_distiller(palette = "Spectral") +
theme_bw()
Next we use the count data to make a map of distinct count sites and save the coordinates of the sites, unique and across all years, for later spatial modeling and plotting.
# make a set of distinct study sites for mapping
site_map <- count_dat %>%
select(circle, easting, northing) %>%
distinct() %>%
select(circle, easting, northing)
Computing a continuous-space model with R-INLA using the
SPDE approach requires construction of four distinct sets of data and
model objects (Blangiardo and Cameletti 2015, Krainski et al. 2018).
First, we create a modeling mesh, which is used to
provide a piecewise linear representation of the continuous spatial
surface, based on a triangulation of the modeled region. Here, the same
mesh will get reused for each of the spatial terms in the model. Second,
we construct an SPDE model object that specifies
properties of the spatial model. Again, we will use the same SPDE object
for each of the spatial terms in the model.
In plain R-INLA, we would also need to create index
vectors and projector matrices (or A matrices as they
are often called). However, in inlabru, these objects are
created automatically, and the user does not need to deal with them
directly. This works by an automatic creation of a
bru_mapper object that knows how to map between mesh nodes
and spatial data locations. Thus, the only two R-INLA
functions needed in the user code of the SPDE modelling steps are
fm_mesh_2d_inla() and inla.spde2.pcmatern(),
and the inla.spde.make.index() and
inla.spde.make.A() functions are only called internally by
the inlabru code itself, or rather
fm_evaluator() is used instead.
There are various things to consider when constructing a mesh (Lindgren and Rue 2015, Blangiardo and Cameletti 2015, Krainski et al. 2018, Bakka et al. 2018). In constructing one, we balance a trade-off between capturing fine-scaled features of the Gaussian random field and computing times. Here, we create two non-convex hulls around the count sites, and then build a triangular mesh by specifying minimum and maximum edge lengths within the inner hull, and within the slightly larger outer hull.
# make a two extension hulls and mesh for spatial model
hull <- fm_extensions(
count_dat,
convex = c(200, 500),
concave = c(350, 500)
)
mesh <- fm_mesh_2d_inla(
boundary = hull, max.edge = c(100, 600), # km inside and outside
cutoff = 50, offset = c(100, 300),
crs = fm_crs(count_dat)
) # cutoff is min edge
# plot it
ggplot() +
gg(data = mesh) +
geom_sf(data = site_map, col = "darkgreen", size = 1) +
geom_sf(data = states, fill = NA) +
theme_bw() +
labs(x = "", y = "")
Next we create an SPDE object to define the model smoothness, with
prior distributions for the variance and range parameters, and the mesh.
Here we assume a Gaussian random field characterized with a Matérn
covariance function with penalized complexity priors (Simpson et
al. 2017) for the practical range (distance where spatial correlation
approaches 0.1) and variation explained by the function (Fuglstad et
al. 2019). The prior for the spatial range is set such that the
probability of a range exceeding 500 km is 0.5. The prior for the
variance explained by the spatial effect is set such that the
probability of a standard deviation exceeding 1 is 0.5 (Krainski et
al. 2018). If one wants to constrain this kind of spatial effect to
integrate to zero, constr=TRUE should be added at this
stage.
# make spde
spde <- inla.spde2.pcmatern(
mesh = mesh,
prior.range = c(500, 0.5),
prior.sigma = c(1, 0.5)
)
In ordinary R-INLA, it would be necessary to construct
projector matrices for each model component, that would include the
different covariate weightings for each spatial model component. In
inlabru, this can instead be handled by either explicitly
multiplying the spatial random fields with the spatial covariates in the
model formula expression, or by specifying the weights
argument for each model component specification. We will here use the
latter approach, since it is the easiest.
In our model, \(\epsilon_s\) is a
spatially varying effect of the log of count effort. Similarly, \(\tau_s\) is a spatially varying effect of
year, so standardized year (with 1987 = 0) is also specified in the
weights argument. \(\alpha_s\) is also an SVC, but it is a
spatially varying intercept. For intercepts it is not necessary to
specify a constant weight of 1. The model component \(\kappa_{s}\) is linked to each observation
site and is not modeled on the mesh.
Note that the current use of the term ‘weights’ is different from
that often encountered when defining mixed effect models in
R. Here it is used to define covariate value
multiplication, as opposed to importance values for
likelihoods in other contexts.
In plain R-INLA, we would need to bundle all the model
data and component projector matrix information using the
inla.stack() function. In inlabru, this is
done automatically, so we skip that step.
The last required input for the analysis is the model formula, which includes information on the prior for explained variation for the unstructured random intercept. We define the prior for \(\kappa_s\) as a penalized complexity prior (Simpson et al. 2017), set such that the probability of the standard deviation associated with the random effect exceeding 1 is 0.01.
# iid prior
pc_prec <- list(prior = "pcprec", param = c(1, 0.1))
Notice that if one wants to constrain a spatial spde effect to
integrate to zero, it should be added constr=TRUE in the
SPDE model definition rather than in the component definition. As we
want to constrain \(\kappa_s\) and it
is a non-spatial term we can use
constr=TRUE in its corresponding label()
component definition below, which imposes a sum-to-zero constraint.
The model described above is translated to inlabru
modeling syntax as:
# components
svc_components <- ~ -1 +
kappa(site_idx, model = "iid", constr = TRUE, hyper = list(prec = pc_prec)) +
alpha(geometry, model = spde) +
eps(geometry, weights = log_hrs, model = spde) +
tau(geometry, weights = std_yr, model = spde)
# formula, with "." meaning "add all the model components":
svc_formula <- count ~ .
Specifying the sf column geometry as input
causes inlabru to extract the spatial coordinate
information for the observation locations from. Alternatively, we could
specify the function st_coordinates instead, but that would
extract the raw coordinates, potentially losing important coordinate
reference system information. Alternatively, we can specify that
CRS-free information explicitly as input, with
cbind(easting, northing) or
st_coordinates(.data.).
Here, we define the response as count, remove the
automatic global intercept with a -1, and then specify the
other terms in the model with label() statements. The first
kappa() statement defines \(\kappa_s\), the site effect, as a normally
distributed and independent (model="iid"), globally
zero-centred (sum to zero constr=TRUE), deviation from
\(\alpha_s\). The second
alpha() statement defines \(\alpha_s\) as a spatially varying intercept
with spatial structure described by the SPDE object called ‘spde’.
Remember that if one wants to constrain this kind of spatial effect to
integrate to zero, constr=TRUE should be added in the SPDE
model definition rather than in the label() arguments. The
third eps() statement defines \(\epsilon_s\) as an SVC for the effect of
count effort, with spatial structure also described by the SPDE object.
The weights for this spatially structured random slope are specified
with the weights=log_hrs argument. The fourth
tau() statement defines \(\tau_s\) as an SVC for the year effect,
with spatial structure described in the SPDE object. The weights for
this spatially structured random slope are specified with the
weights=std_yr argument.
We estimate the model with a call to bru(). First we set
the option to use the (new) experimental way to do internal
computations, see Van Niekerk et. al. (2022), for the sake of computing
speed and better numerics. In the call to bru(), we give
the model components, and specify an observation likelihood with the
model formula for the negative binomial distribution for the counts, and
define the estimation data. Then we ask inla() to compute
WAIC and CPO (temporarily disabled) to evaluate model fit (also, to save
the information necessary for posterior sampling we need
config=TRUE, but this is automatically set by
bru(), as inlabru::predict() relies on it for
posterior sampling and prediction). For computing speed, we choose to
use the simplified Laplace approximation strategy and Empirical Bayes
estimation. The inla() run for this model takes about 1
minutes on a standard laptop computer. Another option is to use the
Variational Bayes approximation as detailed in Van Niekerk and Rue
(2021) and Van Niekerk et. al. (2022).
res <- bru(
svc_components,
like(
formula = svc_formula,
family = "nbinomial",
data = count_dat
),
options = list(
control.compute = list(waic = TRUE, cpo = FALSE),
control.inla = list(int.strategy = "eb"),
verbose = FALSE
)
)
Once computation is complete, we look at the initial results to see how things went. First we check the posterior means for the hyperparameters of the model, mainly the variance components and the spatial ranges of the spatially structured parameters.
# view results
res$summary.hyperpar[-1, c(1, 2)]
## mean sd
## Precision for kappa 2.277815e+00 0.4403796
## Range for alpha 9.530982e+02 271.1085872
## Stdev for alpha 1.934869e+00 0.3719044
## Range for eps 5.203927e+03 2619.2752419
## Stdev for eps 4.150412e-01 0.2266386
## Range for tau 7.424839e+02 244.1578455
## Stdev for tau 6.404605e-02 0.0123772
Next we examine some summaries of the random effect estimates, starting with \(\exp(\alpha_s)\), which is effort-corrected relative abundance at year = 0 (1987), given 1 hour of count effort (i.e., log[1]=0).
summary(exp(res$summary.random$alp$"0.5quant")) # exp(alpha) posterior median
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05803 1.19659 3.45625 6.72359 8.72078 47.49322
Note that to avoid issues due to \(E[\exp(x)|y]\neq\exp[E(x|y)]\), we use the posterior median instead of the posterior mean.
The summary for \(\epsilon_s\) shows variation in the exponent for the effort correction function across space.
summary(res$summary.random$eps$mean) # epsilon
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7731 0.8907 0.9781 0.9669 1.0444 1.1383
The \(\tau_s\) summary shows how long-term, log-linear trends of robin relative abundance have varied across space, from annual decreases of around 10% to annual increases of around 10%.
summary((exp(res$summary.random$tau$"0.5quant") - 1) * 100) # (exp(tau)-1)*100
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -12.936 -5.253 -1.819 -1.220 2.520 11.322
Next we create maps of \(\alpha_s\), \(\epsilon_s\), and \(\tau_s\) to inspect the spatial structure of these parameter estimates. We start by creating a 25-km mapping grid, and then projecting this mapping grid to the modeling mesh.
# get easting and northing limits
bbox <- fm_bbox(hull[[1]])
grd_dims <- round(c(x = diff(bbox[[1]]), y = diff(bbox[[2]])) / 25)
# make mesh projector to get model summaries from the mesh to the mapping grid
mesh_proj <- fm_evaluator(
mesh,
xlim = bbox[[1]], ylim = bbox[[2]], dims = grd_dims
)
Then we populate the mapping grids with parameter estimates (posterior median and range95), turn them into a raster stack, and mask the raster stack to the study area.
# pull data
kappa <- data.frame(
median = exp(res$summary.random$kappa$"0.5quant"),
range95 = exp(res$summary.random$kappa$"0.975quant") -
exp(res$summary.random$kappa$"0.025quant")
)
alph <- data.frame(
median = exp(res$summary.random$alpha$"0.5quant"),
range95 = exp(res$summary.random$alpha$"0.975quant") -
exp(res$summary.random$alpha$"0.025quant")
)
epsi <- data.frame(
median = res$summary.random$eps$"0.5quant",
range95 = (res$summary.random$eps$"0.975quant" -
res$summary.random$eps$"0.025quant")
)
taus <- data.frame(
median = (exp(res$summary.random$tau$"0.5quant") - 1) * 100,
range95 = (exp(res$summary.random$tau$"0.975quant") -
exp(res$summary.random$tau$"0.025quant")) * 100
)
# loop to get estimates on a mapping grid
pred_grids <- lapply(
list(alpha = alph, epsilon = epsi, tau = taus),
function(x) as.matrix(fm_evaluate(mesh_proj, x))
)
# make a terra raster stack with the posterior median and range95
out_stk <- rast()
for (j in 1:3) {
mean_j <- cbind(expand.grid(x = mesh_proj$x, y = mesh_proj$y),
Z = c(matrix(pred_grids[[j]][, 1], grd_dims[1]))
)
mean_j <- rast(mean_j, crs = epsg6703km)
range95_j <- cbind(expand.grid(X = mesh_proj$x, Y = mesh_proj$y),
Z = c(matrix(pred_grids[[j]][, 2], grd_dims[1]))
)
range95_j <- rast(range95_j, crs = epsg6703km)
out_j <- c(mean_j, range95_j)
terra::add(out_stk) <- out_j
}
names(out_stk) <- c(
"alpha_median", "alpha_range95", "epsilon_median",
"epsilon_range95", "tau_median", "tau_range95"
)
out_stk <- terra::mask(
out_stk,
terra::vect(sf::st_union(states)),
updatevalue = NA,
touches = FALSE
)
Finally, we plot the SVCs with the following code. We plot the posterior median and 95% uncertainty width (“range95”) for \(\exp(\kappa_s)\), \(\exp(\alpha_s)\), \(\epsilon_s\), and \(100(\exp(\tau_s)-1)\).
make_plot_field <- function(data_stk, scale_label) {
ggplot(states) +
geom_sf(fill = NA) +
coord_sf(datum = NA) +
geom_spatraster(data = data_stk) +
labs(x = "", y = "") +
scale_fill_distiller(scale_label,
palette = "Spectral",
na.value = "transparent"
) +
theme_bw() +
geom_sf(fill = NA)
}
make_plot_site <- function(data, scale_label) {
ggplot(states) +
geom_sf() +
coord_sf(datum = NA) +
geom_sf(data = data, size = 1, mapping = aes(colour = value)) +
scale_colour_distiller(scale_label, palette = "Spectral") +
labs(x = "", y = "") +
theme_bw() +
geom_sf(fill = NA)
}
# medians
# fields alpha_s, epsilon_s, tau_s
pa <- make_plot_field(
data_stk = out_stk[["alpha_median"]],
scale_label = "posterior\nmedian\nexp(alpha_s)"
)
pe <- make_plot_field(
data_stk = out_stk[["epsilon_median"]],
scale_label = "posterior\nmedian\nepsilon_s"
)
pt <- make_plot_field(
data_stk = out_stk[["tau_median"]],
scale_label = "posterior\nmedian\n100(exp(tau_s)-1)"
)
# sites kappa_s
ps <- make_plot_site(
data = cbind(site_map, data.frame(value = kappa$median)),
scale_label = "posterior\nmedian\nexp(kappa_s)"
)
# range95
# fields alpha_s, epsilon_s, tau_s
pa_range95 <- make_plot_field(
data_stk = out_stk[["alpha_range95"]],
scale_label = "posterior\nrange95\nexp(alpha_s)"
)
pe_range95 <- make_plot_field(
data_stk = out_stk[["epsilon_range95"]],
scale_label = "posterior\nrange95\nepsilon_s"
)
pt_range95 <- make_plot_field(
data_stk = out_stk[["tau_range95"]],
scale_label = "posterior\nrange95\n100(exp(tau_s)-1)"
)
# sites kappa_s
ps_range95 <- make_plot_site(
data = cbind(site_map, data.frame(value = kappa$range95)),
scale_label = "posterior\nrange95\nexp(kappa_s)"
)
# plot together
multiplot(ps, pa, pe, pt, cols = 2)
The map for the posterior mean of \(\tau_s\) shows that robins have decreased
in the southern part of the study area and increased in the northern
part. This demonstrates how the wintering range of robins is shifting
northward as winters become warmer due to climate change.
# plot together
multiplot(ps_range95, pa_range95, pe_range95, pt_range95, cols = 2)